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GEOSTATISTICAL ANALYSIS OF RESOURCE ISLANDS UNDER 
ARTEMISIA TRIDENTATA IN THE SHRUB-STEPPE 


Jonathan J. Halvorson!, Harvey Bolton, Jr.2, Jeffrey L. Smith®?, and Richard E. Rossi? 


ABSTRACT.—Desert plants ean inflttenee the pattern of resources in soil resulting in small-scale enriched zones. 
Although conceptually simple, the shape, size, and orientation of these “resource islands” are difficult to study in detail 
using conventional sampling regimes. To demonstrate an alternative approach, we sampled soil under and around indi- 
vidual Artemisia tridentata (sagebrush), a dominant shrub of cool desert environments, and analyzed the samples with 
univariate statistics and geostatisties. Univariate statisties revealed that soil variables like total inorganie-N, soluble-C, 
and mierobial biomass-C were distributed with highest mean values within about 25 cm of the plant axis and significant- 
ly lower mean values at distances beyond 60 em. However, such simple analyses restricted our view of resource islands 
to identically sized, symmetrical accumulations of soil resourees under each plant. 

Geostatisties provided additional information about spatial characteristies of soil variables. Variography revealed that 
samples separated by a distance of less than about 70 em were correlated spatially. Over 75% of the sample variance was 
attributable to spatial variability. We modeled these spatial relationships and used kriging to predict valties for unsam- 
pled locations. Resulting maps indicated that magnitude, size, and spatial distribution of soil resource islands vary 
between individual plants and for different soil properties. Maps, together with cross-variography, further indicate that 
resouree islands under A. tridentata are not always distinguishable from the surrounding soil by sharp transition bound- 
aries and may be asymmetrically distributed around the plant axis. 


Key words: resource islands, geostatistics, Artemisia tridentata, nutrient availability, kriging, spatial correlation. 


Recognition that individual plants can sig- 
nificantly affect the local soil environment 
dates back to at least the middle of the nine- 
teenth century (see Zinke 1962) and has been 
documented for many plant forms including 
broadleaf and coniferous trees (Zinke 1962, 
Everett et al. 1986, Doescher et al. 1987, 
Belsky et al. 1989), bunch grasses (Hook et al. 
1991), herbaceous legumes (Halvorson et al. 
1991), and, in particular, desert shrubs (e.g., 
Fireman and Hayward 1952, Garcia-Moya and 
Mckell 1970, Nishita and Haug 1973, Barth and 
Klemmedson 1978, Burke 1989, Burke et al. 
1989, Virginia and Jarell 1983, Bolton et al. 
1990, 1993). Soil associated with plants typi- 
cally contains greater concentrations of limit- 
ing resources (e.g., N, P), contains larger pop- 
ulations of soil microorganisms, and exhibits 
higher rates of nutrient cycling processes like 
mineralization (Charley and West 1977, Bolton 
et al. 1990) and denitrification (Virginia et al. 
1982). These small-scale enriched zones, vari- 
ously termed “fertile islands” (Garner and Stein- 
berger 1989), “isles of fertility” (West 1981, 


Whitford 1986), “resource islands’ (Reynolds 
et al. 1990), or “ecotessara’ (Jenny 1980), are 
hypothesized to result from several mecha- 
nisms (Garner and Steinberger 1989), includ- 
ing litter-fall or stemflow (Zinke 1962), 
decreased erosion or increased deposition 
(Coppinger et al. 1991), microclimatological 
amelioration of the soil (Smith et al. 1987, 
Pierson and Wight 1991), or inputs of resources 
via insects, birds, or animals (Davidson and 
Morton 198-4). 

Detailed knowledge of the size and internal 
dynamics of resource islands is important for 
understanding energy flux, mass transport, 
and nutrient cycling processes at the scale of 
the individual plant. Resource islands may also 
connote a tier in a progressive, hierarchical 
mosaic of plant and animal habitats, resource 
distributions, and biogeochemical processes 
(i.e., patches sensu Kotliar and Wiens 1990). 
Estimates of the distribution and numbers of 
resource islands in the landscape may aid in 
understanding population level processes and 
can be used to refine regional estimates of 
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energy flow and mass transfer. Furthermore, 
interrelationships of large numbers of individ- 
ual resource islands may influence ecosystem 
structure, function, and stability (Reynolds et al. 
1990, Schlesinger et al. 1990, Halvorson et al. 
Lot 

Although conceptually simple, the size, 
shape, and orientation of resource islands are 
not easy to evaluate. Previous studies have 
typically been based on relatively small num- 
bers of samples collected using a binary regime 
(i.e., samples collected beneath the plant versus 
samples collected “away” from the plant) or 
along a transect passing from plant to bare soil. 
Such an approach cannot be used to provide a 
detailed spatial analysis of resource concentra- 
tions or processes in the soil that are likely to 
exhibit complex responses to landscape and 
microsite variations (Burke et al. 1989). 
Additionally, data collected from different 
locations (or depths) have often been analyzed 
using inferential statistics such as ANOVA or 
t tests that assume samples are spatially inde- 
pendent and identically distributed. However, 
these assumptions may be dubious, if untested, 
since ecological phenomena are often spatially 
or temporally correlated and their frequeney 
distributions are rarely normal (Rossi et al. 
1992). Recently, a branch of applied statistics, 
known as geostatistics, has been demonstrated 
to be useful for determining spatial correlations 
among ecological data and for estimating values 
at unsampled locations (e.g., Robertson 1987, 
Robertson et al. 1988, Rossi et al. 1992). 

Objectives of this study were (1) to use geo- 
statistics to describe and model the spatial 
continuity of soil variables around individual 
plants, (2) to use this information to produce 
graphical representations or maps of specific 
resource islands, and, finally, (3) to quantify 
the spatial correlation between plants and soil 
variables. We examined Artemisia tridentata 
(sagebrush), a prominent shrnb of cool desert 
environments (West 1983) previously known 
to affect the distribution of resources in the 
soil. Several workers have measured higher 
coneentrations of resourees such as total-C, 
total-N, inorganic-N, and higher rates of N 
evcling in soil beneath A. tridentata than in 
nearby open soil using a binary sampling 
regime (e.g., Burke 1989, Burke et al. 1989, 
Bolton et al. 1990, 1993). However, these stud- 
ies have not accounted for possible spatial auto- 
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correlation of the samples, evaluated resource 
islands of individual plants, nor quantified the 
scale of soil heterogeneity beneath A. tridentata. 

Geostatistics has previously been used to 
describe environmental and soi] parameters 
associated with A. tridentata. For example, 
Pierson and Wight (1991) used one-dimen- 
sional geostatistics to analyze spatial and tem- 
poral variability of soil temperature under A. 
tridentata. Halvorson et al. (1992) demonstrat- 
ed that geostatistics was an appropriate means 
of measuring resource islands at the scale of 
an individual A. tridentata plant. Jackson and 
Caldwell (1993a) attempted to quantify the 
scale of nutrient heterogeneity around indi- 
vidual A. tridentata and Pseudoroegneria spi- 
cata in a native sagebrush steppe using semi- 
variograms. They demonstrated increasing 
autocorrelation of soil nutrients at spatial 
scales <1 m but did not determine whether 
small-scale effects were attributable to individ- 
ual plants or an artifact of the nested sampling 
design used. More recently, they constructed 
kriged maps that showed relatively high con- 
centrations of soil variables like soil organic 
matter, extractable phosphate, and potassium 
near Pseudoroegneria tussocks but not Artemisia 
shrubs (Jackson and Caldwell 1993b). However, 
these kriged maps did not directly quantify 
spatial covariation between locations of indi- 
vidual plants and resource islands. Further, 
Jackson and Caldwell did not observe autocor- 
relation for microbial processes at any scale 
that was measured. 

To meet our objectives, we applied geosta- 
tistics in three steps. First, we characterized 
and modeled the similarity between samples 
as a function of their separation distance and 
direction. Second, we used this relationship to 
interpolate values at unsampled locations 
directly under and near individual plants. 
Finally, we quantified spatial covariation be- 
tween soil properties and plants. 


STUDY SITE 


The study was conducted at the Arid Land 
Ecology (ALE) Reserve, located on the 
Hanford Site in south central Washington (see 
Bolton et al. 1990 for details). There, remnants 
of the native Artemisia tridentata—Elytrigia 
spicata association occur on silt loams of the 
Warden or Ritzville series. This perennial 
shrub-steppe is the largest grassland-type in 
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North America and covers more than 640,000 
km? of the Intermountain Pacific Northwest 
too dry to support forests (Daubenmire 1970, 
Rogers and Rickard 1988). {1 an undisturbed 
state the A. tridentata—E. spicata association 


would be composed typically of three layers of 


vegetation: an overstory shrub (Artemisia tri- 
dentata tridentata), a large caespitose perenni- 
al grass (Elytrigia spicata [formerly Agropyron 
spicatum|), and a small cacspitose perennial 
grass (Poa secunda) growing on soil with a thin 
cryptogamic crust (Daubenmire 1970). 
However, following disturbance such as 
tillage, grazing, or fire, the alien annual grass 
Bromus tectorum becomes established. 


METHODS 


Soil Collection and Analysis 


Cores of surface soil (10.5 em dia. X 5 cm 
deep) were collected at 41 specific locations 
within 2 X 2-m plots centered on mature A. 
tridentata individuals (Fig. 1). Samples were 
located so as to minimize the number of data 
points needed for analysis of spatial character- 
istics and to avoid preferential clustering. We 
sampled five identically oriented plots (205 
points) in March 1991, when levels of soil 
moisture and microbial biomass activity were 


North (m) 





Fig. 1. Schematic of a typical sampling plot. Each plot 
(five total) was centered on an Artemisia tridentata plant. 
Dashed circles show the location of 41 soil cores (10.5 cm 
dia. X 5 cm deep). All five plots were oriented as shown. 
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high. All plots were located within approxi- 
mately 20 m of each other within a flat area 
with randomly spaced plants. Multiple plots 
were sampled for two reasons: first, to assess 
spatial characteristics of resource islands by 
basing our calculations on several examples 
rather than a single instance; and, second, to 
provide replicates in the event that no spatial 
dependence of soil properties was observed. 
Data from all plots were combined to consider 
spatial dependence of soi] properties around 
several A. tridentata plants simultaneously. 
This approach was chosen because it provided 
a more generalized evaluation of resource 
islands under individual A. tridentata and 
greatly increased the number of data pairs at 
anv separation clistance. 

Estimates of plant location were required 
for cross variography (see below). Therefore, 
vegetation maps were produced from vertical 
photographs. Each plot was divided into 5 x 
5-cm squares. Each square was classified into 
one of three groupings—bare, grass species, 
or A. tridentata—based on predominant cov- 
erage. For this work no attempt was made to 
distinguish among grass species. 

Each soil sample was sieved (5 mm), mixed, 
and analyzed for a variety of soil variables. For 
this work we present data only for water-solu- 
ble forms of C, total inorganic-N (i.e., nitrate 
+ ammonium), and soil microbial biomass-C. 
Soluble soil C (Hj0-C) was extracted with 
room temperature deionized water and ana- 
lyzed using an infrared gas analyzer (Ionics 
Inc., Watertown, Massachusetts). Total inor- 
ganic nitrogen (TI-N) was extracted within 48 h 
of collection from 10-g subsamples of soil 
using 25 ml 2M KCl and analyzed colorimetri- 
cally (Alpkem Corp., Clackamas, Oregon). Soil 
microbial biomass-C (SIR-C) was estimated 
from the respiratory response of soil to glucose, 
a source of C readily utilized by heterotrophic 
soil microorganisms (Anderson and Domsch 
1978). Ten-gram samples of soil were placed 
in 40-ml glass vials, moistened with deionized 
H30, covered with Parafilm, and incubated in 
the dark at 23.5 + 0.5°C for I wk. Each sam- 
ple was then amended with a glucose solution 
at the rate of 600 mg glucose (240 mg C) kg} 
soil, bringing the final H,O content of the soil 
to 20-25% (w/w; equivalent to 30-50 kPa). 
Glass vials were sealed with silicone septa and 
incubated for 3 h. Soil respiration was mea- 
sured by gas chromatography and related to 
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estimates of soil microbial biomass-C with 
equations developed by Anderson and Domsch 
(1978). 


Univariate Statistics 


Univariate statistics were calculated for 
each soil parameter. For classical inferential 
statistics, data were also assigned to one of five 
distance classes depending upon sample loca- 
tion within a plot. These classes can be envi- 
sioned as concentric rings located at increas- 
ing distances from the center of the plot. The 
first distance class comprised samples collected 
directly under A. tridentata (distance = 0 cm, 
n = 1 per plot), followed by the second 
(approximate distance = 25 cm, n = 8 per 
plot), third (approximate distance = 60 cm, n 
= § per plot), fourth (approximate distance = 
110 cm, n = 12 per plot), and fifth (approximate 
distance = 130 cm, n = 12 per plot). Average 
values of soil properties in each distance class 
were plotted as a function of radial distance 
from the plant axis. Following variography (see 
below), logy) transformed samples deemed 
spatially independent were compared with 
ANOVA using plot as a blocking factor. 


Geostatistics 


VARIOGRAPHY.—We evaluated spatial char- 
acteristics of each soil parameter with the non- 
ergodic autocorrelation function (Srivastava 
and Parker 1989) and summarized results 
graphically as correlograms. Like variograms, 
correlograms represent the average degree of 
similarity hetween samples as a function of 
their separation distance (lag) and direction. 
Unlike the variogram, the correlogram filters 
out the effects of changes in both lag means 
and lag variances. Each point in a correlogram 
was calculated from this equation: 


1 
p SN. 


Nth) 
i=11[z(x;) — m_),|[z(x; +h) — m4), 


S_h eet 


where z(x;) and z(x; + h) are two data points 
separated by the distance (lag) h. Datum z(x;) 
is the tail and z(x; + h) is the head of the vec- 
tor, N(h) is the total number of data pairs sepa- 
rated by lag h, m_, and m,}, are means of the 
points that correspond to tail and head of the 
lag, respectively, and S_;, and S,), are standard 
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deviations of tail and head values of the lag, 
respectively. We chose the correlogram 
because it removes the effects of lag means 
and standardizes by the lag variances (Rossi et 
al. 1992). For this work we express correlo- 
grams in the form of a standardized variogram 
by subtracting each p*(h) from 1 (saaks and 
Srivastava 1989, Rossi et al. 1992). 

Correlograms were first calculated solely as 
a function of lag distance (i.e., the omnidirec- 
tional case) without considering any differ- 
ences in spatial continuity with direction (i.e., 
anisotropy). However, since resource islands 
need not be symmetric (e.g., Zinke 1962), we 
also ealculated directional correlograms. For 
each soil property a separate correlogram was 
calculated for samples oriented 0°, 45°, 90°, 
and 135° (£15° tolerance) from each other. 
Since correlograms are symmetric about the 
origin (.e,, 0° = 180°) 415° = 225 7 cto 
45°, 90°, and 135° directions correspond to 
samples aligned along east-west, northeast- 
southwest, north-south, and northwest-south- 
Cast anes: 

Because the data of each plot were concate- 
nated during this analysis, local anisotropies 
(i.e., anisotropies specific to each plot) were in 
effect combined. Thus, any directional effeets 
we observed were a composite of the five plots 
and presumably indicative of overall direction- 
al trends. To identify directions of maximum 
and minimum continuity, we estimated the lag 
distance corresponding to a common value for 
each directional correlogram (Isaaks and 
Srivastava 1989). The directional correlogram 
with the greatest lag associated with a correlo- 
gram value of 1 was identified as the direction 
of greatest continuity. The correlogram with the 
smallest lag corresponding to 1 was deemed 
the dircetion of minimum continuity. 

The empirically determined scatter of data 
points in each correlogram was fit with models 
known to produce a positive definite kriging 
system (i.e., matrices that provide both a unique 
solution and a positive estimation variance; 
Isaaks and Srivastava 1989). Such models typi- 
cally contain several salient features known as 
nugget, range, and sill. The nugget is the 
amount of variance not explained or modeled 
as spatial correlation. It is the apparent ordi- 
nate intercept and is due to (1) unsampled cor- 
relation below the smallest lag and (2) experi- 
mental error (Rossi et al. 1992). A small nugget 
relative to the sill indicates that a large pro- 
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portion of the sample variability is modeled as 
spatial dependence. Conversely, a large nugget 
indicates less sample variability can be mod- 
eled as spatial dependence. The sill is charac- 
terized by a leveling off of the correlogram 
model. If present, it indicates that spatial cor- 
relation is, on average, constant. However, if 
spatial correlation continues to change at lags 
greater than those considered in the correlo- 
gram, then a sill will not be apparent. The lag 
value when the correlogram model reaches 
the sill is known as the range. It represents 
the maximum separation distance within 
which samples are spatially correlated. At lags 
> the range, the sill of the variogram may 
approach the sample variance (Barnes 1991). 

KRIGING TO ESTIMATE DATA AT UNSAMPLED 
LOCATIONS.—Kriging has been likened to 
“multiple linear regression with a few twists” 
(Rossi 1989). In classical multiple linear regres- 
sion, an estimate of the dependent variable, Y, 
is calculated from a weighted linear combina- 
tion of independent variables where each is 
measured at about the same location in time 
or space. Usually, only a single value of Y is 
estimated. Similarly, in kriging, z*(x,), the 
estimated value of the variable for an unsam- 
pled location (x,), is calculated as a weighted 
linear combination of the surrounding sam- 
pled neighbors: 


N 
a a a(x) (2) 


where the z(x;)'s are the sampled values at 
their respective locations, and the g;s are the 
weights associated with each sample value. In 
ordinary kriging, weights used to estimate 
z*(x,) are chosen so that the resulting estimate 
is unbiased and has a minimum estimation 
variance and sum to unity. kriging incorpo- 
rates a model of spatial continuity (here the cor- 
relogram model) and accounts for the degree 
of clustering of nearby samples and their dis- 
tance to the point being estimated (Isaaks and 
Srivastava 1989). We used ordinary point kriging 
to estimate values of soil properties at unsam- 
pled locations. For each plot we estimated val- 
ues for the nodes of a 5 X 5-cm?2 grid. Each 
predicted value was based on a minimum of 6 
and a maximum of 12 neighbors located with- 
in a 0.8-m circular search radius. 
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CROSS-VARIOGRAPHY.—In addition to spatial 
characteristics of single soil properties, we 
also determined how soil properties covaried 
with plants. We modeled spatial covariation 
with p*,p(h), the estimated nonergodic cross- 
correlogram. Like the correlogram, it accounts 
for both variables’ fluctuating lag means and 
variances (Isaaks and Srivastava 1989, Rossi et 
al. 1992). Because comparisons between contin- 
uous variables (i.e., TI-N, SIR-C, and H,O-C 
data ) and discrete variables (i.e., plant data) 
might be complicated by a “contact effect” 
(Luster 1985), we converted TI-N, SIR-C, and 
H,O-C data to binary variables using an indi- 
cator transformation (Journel 1983). For this 
work, continuous data values of TI-N, SIR-C, 
and H,O-C were coded 1 if they were greater 
than the local (within-plot) median, or 0 follow- 
ing Halvorson et al. (in review). Cross-correlo- 
grams were then calculated for grass species 
or A. tridentata and indicator transformed TI-N, 
SIR-C, and H,0-C data using the equation, 





le (x;,Z4) — ma_, tl%.28) - mp, (3) 


S448 Bh 


where N(h) is the total number of data pairs 
separated by vector h, L4(x;,z,4) is the coded 
plant data (equal to 1 if the specified plant 
type was present or 0 if absent) at some loca- 
tion (x;), ma_, and Sq , are the mean and 
standard deviation, respectively, for the plant 
variable at those data locations that are —h 
away from a soil property data location. 
Similarly, Ip(x,,zp) is the coded soil variable 
data (equal to 1 if the data value is greater 
than the local plot median or else 0) at location 
(x},), Mp, and Sp, are the mean and standard 
deviation of the soil variable indicator calcu- 
lated for those locations that are +h away 
from a plant variable data location. Note, 
when h is 0, equation 3 is equivalent to the 
Pearson correlation coefficient (Isaaks and 
Srivastava 1989). 

Unlike the correlogram, values calculated for 
the cross-correlogram may not be symmetric 
about the origin because both the order and 
direction are switched when variables are re- 
versed (Isaaks and Srivastava 1989). Conse- 
quently, we calculated individual cross-correl- 
ograms for the 0°, 45°, 90°, 135°, 180°, 225°, 
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270°, and 315° directions (+30° tolerance). 
These correspond to soil samples aligned to 
the east, northeast, north, northwest, west, 
southwest, south, or southeast of a plant. 


RESULIrS 


Summary statistics indicated that samples 
of TI-N and H,0-C were positively skewed, 
while samples of SIR-C were more normally 
distributed (Figs. 2A—C). Total inorganic-N 
ranged from 0.6 mg / kg soil to a maximum of 
23.6 mg / ke soil (Fig. 2A). The mean value for 
TI-N of 3.8 mg / kg soil compared reasonably 
to the values reported by Bolton et al. (1990) of 
4.1 and 4.9 mg / kg soil for open soil crust and 
A. tridentata soil, respectively. Values observed 
for H,O-C ranged widely from 9.8 mg / kg soil 
to 633.9 mg / kg soil (Fig. 2B). Estimates of 
SIR-C ranged from less than 200 to over [SOO 
mg / kg soil (Fig. 2C). The average value for 
SIR-C, 750 mg/ kg, was within the range 
reported by Burke et al. (1989) and equivalent 
to about 980 kg C / ha soil assuming a bulk 
density of 1.3 (Bolton et al. 1990). Comparative- 
ly, Smith and Paul (1990) reported average 
microbial biomass pool size for grassland sys- 
tems of 1090 kg C /ha. 

Univariate statistics also indicated how soil 
properties varied with distance from the A. 
tridentata axis (i.c., center of the plot; Figs. 
3A-C). Concentrations of HygO-C and SIR-C 
were greatest within 25 cm of the plant axis 
and lowest at distances bevond 60 cm (Figs. 
3B, C). A similar pattern was observed for TI-N 
except that mean coneentration was low in soil 
collected from directly beneath the A. trideu- 
tata plant (Fig. 3A) and from distanees beyond 
60 cm. This somewhat unexpected finding of a 
resouree “hole” in the center of the resource 
island may be indieative of differences in the 
cycling of N under sagebrush and grass plants. 

Variography indicated that samples of TI-N, 
H,O-C, and STR-C were spatially correlated 
(Fig. 4). Correlograms for SIR-C and TI-N 
exhibited similar ranges of about 0.7 or 0.8 m. 
The correlogram for HyO-C was similar to the 
others at small lag distances and equaled the 
sample variance at a range near 0.7 m. How- 
ever, at greater lags, correlogram values for 
H.O-C increased above the sample variance 
and did not appear to reach a sill until lags 
were greater than 1 m. A correlogram sill 
greater than I for H»O-C can oecur if the 
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Fig. 2. Frequency histograms and some basic summary 
statistics for (A) total inorganic-N (TI-N), (B) water solu- 
ble-C (H1,0-C), and (C) soil microbial biomass-C (STR-C). 


majority of sample values are collected from 
an area with dimensions equal to or less than 
the variogram range (Barnes 1991) or if dis- 
crete regions of high and low concentration 
occur at lags greater than the maximum lag in 
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Fig. 3. Mean + standard error for (A) TI-N, (B) H,0-C, 
and (C) SIR-C. Numbers in parentheses are the number 
of data points (from five plots combined) that contribute 
to each estimate. 
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Fig. 4+. Omnidirectional correlograms for TI-N, H,0-C, 
and SIR-C. Each point shown was calculated from a mini- 
mum of 253 pairs (range 253-644). 


the correlogram (i.e., an incompletely mod- 
eled “hole” effect). The apparent nuggets for 
all three soil parameters suggested that more 
than 75% (a correlogram value of 0.25 or less) 
of total sample variability could be modeled as 
spatial dependence. 

Ranges observed in correlograms of soil 
properties were used to establish the separation 
distance beyond which correlation between 
samples could not be distinguished from sam- 
ple variance. In other words, samples separat- 
ed by distances greater than the range were 
candidates for analysis using more traditional 
statistical techniques that assume indepen- 
dence such as ANOVA. For our data, samples 
in the first two distance classes were thus 
combined and compared to samples from the 
last two classes because they were separated by 
more than about 0.8 m. Samples in the third 
distance class, lying at an intermediate dis- 
tance between the center and outer boundary 
of the plot, were excluded from analysis. 
Analysis of variance of logjg transformed data 
using the five plots as blocking factors showed 
mean concentrations of TI-N, H,O-C, and 
SIR-C to be significantly greater (P < .001) 
within 29 cm of the plant axis than values col- 
lected = 1.07 m away from the plant axis (Table 

Omnidirectional correlograms character- 
ized spatial correlation or continuity purely as 
a function of lag distance. However, by consid- 
ering the orientation of samples in addition to 
their lag distance, directional anisotropies were 
suggested. Directional correlograms showed 
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TABLE 1. Sumiuary statistics and randomized complete block ANOVAs for logy, transformed T1-N, H,O0-C, and SIR- 
C. Data are summarized into two sample location classes: Near (all measurements collected from within 29 em of the 
plant axis; 2 = 45) and Away (samples collected at distances greater than 107 cm from the plant axis; n = 120). Average 
separation distance between the two location classes was 99.4 em. 


Mean Standard deviation Standard error 
Near Away Near Away Near Away 
Tian O.5S 0.39 0.34 0.32 0.05 0.03 
[ IoC )-C 2.01 1.49 0.31 0.26 0.05 0.02 
SIR-C 2.93 oc 0.14 OF 0.02 0.02 
Source of variation df MS F if 
Plot 4 0.59 6.14 <.001 
TI-N Sample Location I 1.14 11.81 001 
Error 159 0.10 
Plot 4 0.2] 2a 026 
I1,0-C Sample Location ] 8.79 22.49 <.001 
Error 159 0.07 
Plot 4 0.06 ee 060 
SIR-C Sample Location | 0.36 13.63 <.00]1 
Error 159 0.03 


differences in both nuggets and ranges (Fig. 5). 
Generally, the largest apparent nuggets were 
observed in correlograms oriented in the 0° 
(east-west) and 45° (northeast-southwest) 
directions. With the exception of H,O-C, 
these correlograms had estimated nuggets of 
0.4 or more. Conversely, correlograms calcu- 
lated for the 90° (north-south) and 135° 
(northwest-southeast) directions generally 
exhibited nuggets of < 0.2. 

Directions of maximum and minimum con- 
tinuity were identified. For TI-N, maximum 
continuity was observed in the 45° direction 
while lower but similar ranges of continuity 
were observed in the other three directional 
correlograms (Fig. 5, left column). Little 
anisotropy was observed in directional correl- 
ograms for H,O0-C, suggesting that spatial cor- 
relation could be adequately modeled with a 
single isotropic correlogram (Fig. 5, center 
column). Like TI-N, the direction of maximum 
continuity observed for SIR-C was 45°, with a 
direction of minimum continuity in the 135° 
direction (Fig. 5, right column). The aniso- 
tropies for TI-N and SIR-C were accounted 
for in kriging by using a model that evaluated 
both distance and direction (Table 2). 

Maps of the estimates generated with these 
models and ordinary kriging were constructed 
for each soil property in each plot. Taken 
together they suggest that gencralizations 
about spatial distribution of resources in the 
soil beneath A. tridentata can be complicated 


by the variation observed between individual 
plants and specific soil properties. For exam- 
ple, distinct “islands” of TI-N were not alwavs 
clearly associated with A. tridentata. Instead, 
in three of five plots highest coneentrations of 
TI-N appeared to be associated with location 
of grasses (Fig. 6, plots A,B,E). In plots C and 
D, concentrations of TI-N were highest in the 
vicinity of the A. tridentata canopy. However, 
vegetation maps of these plots indicate that 
grass species were present near the A. triden- 
tata plant. In plot E smallest concentrations of 
TI-N were predicted to lie under the A. tri- 
dentata plant. 

Conversely, highest accumulations of H,O- 
C were clearly associated with A. tridentata in 
kriged maps. In each plot high concentrations 
of H,O-C were localized near the plot center 
under the plant canopy. Loeation of grass 
species did not always appear to coincide 
strongly with high concentrations of H ,O-C 
(see NW 1/4 of plot A, NE and SE 1/4 of plot 
B, and NE I/4 of plot E). However, high con- 
centrations of both TI-N, and H,O-C did 
coincide with location of grasses in the NW 
1/4 of plot B and SW 1/4 of plot E. 

Maps of SIR-C indicate that islands of soil 
biomass-C were present under plants but not 
apparently specific to any particular type of 
vegetation. Relatively high concentrations of 
SIR-C were estimated under A. tridentata 
plants in all plots. However, SIR-C was also 
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Fig. 5. Directional correlograms for TI-N, H,O-C, and SIR-C in the 0°, 45°, 90°, and 135° directions calculated with 
a tolerance of 15°. Each point shown was calculated from a minimum of 47 pairs (range 47-212). Vertical lines corre- 
spond to a correlogram value of 1 and were used to identify directions of maximum and minimum spatial continuity. 
Open symbols are estimates of the apparent nugget (the apparent ordinate) fit by eve. 


accumulated elsewhere in relation to the loca- 
tion of grass species. High concentrations of 
SIR-C were observed in several instances not 
associated with high concentrations of either 
TI-N or H,O0-C (e.g., NE and SW 1/4 of plot 
B, NW 1/4 of plot D). 

Cross-variography indicated how TI-N, 
H,O-C, and SIR-C covaried spatially with 
respect to A. tridentata and grass species. 
Indicator transformed TI-N data were similarly 
and positively correlated with A. tridentata and 
grass species at a lag of 0 (equivalent to the 
Pearson correlation coefficient; Figs. 7A,B). 
However, correlation varied with increasing lag 
distance (i.e., showed spatial dependence) in a 
different manner for each. For A. tridentata 
highest positive correlations with TI-N were 


observed in the 45° (northeast) and 90° (north) 
directions (Fig. 7A). The range over which A. 
tridentata remained positively correlated with 
TI-N was longest in the 45° and 90° directions, 
extending to about I m and 0.75 m, respec- 
tively. Positive correlations with TI-N were 
observed for other directions too but only to a 
lag of about 0.5 m. In contrast, grass species 
were positively correlated to above-median 
TI-N concentrations in the 315° (southeast), 
270° (south), and 225° (southwest) directions 
(Fig. 7B). The range over which grass species 
were positively correlated to TI-N in these 
directions was Jess than that for A. tridentata. 
In other directions grass species were uncor- 
related or negatively correlated to TI-N at lags 
above 0. 


aes 


TABLE 2. Model!.? parameters used for ordinary kriging 
ALE fr properties. 


Soil parameter 


TI-N 
l-p*yx = 154 + 0.425 Sph(0.27) + 0.524 Sph(1.80) 
L-p*) a= = 154 + 0.428 Sph(0.71) + 0.524 Sph(1.20) 
STR-C 
l-p*ys = .072 + 0.381 Sph(0.23) + 0.717 Sph(1.87) 
I-p* 135 = .072 + 0.351 Sph(0.74) + 0.717 Sph(0.75) 
11,0-( 
]-p* = [40 Sphi1.3) 


'Nfodels shown for TI-N and SIR-C are a combination of a nugget constant 


anc tvo spherical models. The spherical model, denoted “Sph.” is an autho- 


rized model commonly nsed to model variograms. The number that precedes 
“Sph” can be thought of as the local sill for that model while the number in 
parentheses is the range at which the local sill is reached (see Isaaks and 
Srivastava 1989). For a correlogram, the standardized form is 1l-p*(h) = 1.5 
lag/range) 0.5 (lag/range)? if lag < = range, else 1 if otherwise 

2NoTE: Geostatisticians often distinguish between the nngget used for diag- 
nostic purposes, whieh is the apparent ordinate intercept, and the nugget 
valne, which is used in modeling. 


Indicator transformed H,O-C data were 
positively correlated with A. tridentata but not 
with grass species at a lag of 0 (Figs. 7C,D). As 
with TI-N, highest correlations with A. triden- 
tata were observed in the 45° and 90° direc- 
tions (Fig. 7C). Similar patterns of spatial 
dependence were observed for all directions. 
The distance to which Hl,O-C remained posi- 
tively correlated with A. tridentata ranged 
from a minimum of about 0.5 m to a maximum 
of near 0.75 m in the 45° and 90° directions. 

Unlike A. tridentata, H2O-C was not positive- 
ly correlated with grass species (Fig. 7D). 
Instead, HgO-C was moderately negatively 
correlated in the 0°, 45°, 90°, and 135° diree- 


tions, meaning grass species were more asso- 


ciated with below-median concentrations of 


H5O-C. At lags greater than about 0.2 m, little 
change in cross-correlograms was observed, 
indicating only a weak spatial dependence. 

In contrast to H,O-C, cross-correlograms 
indicated that SIR-C was slightly more corre- 
lated with grass species than with A. tridenta- 
ta at a lag of 0 (Figs. 7E,F). Like other soil 
properties, strongest positive correlation 
between A. tridentata and SIR-C was 
observed in the 45° direction, which also 
remained positively correlated to lags in 
excess of 1 m (Fig. 7E). Lowest correlations 
with A. tridentata were to the 270° and 225 
directions. Indicator transformed SIR-C data 
were most correlated to grass species in the 
225), 270m. andealo dircetioncs bigs 7h). 


aeatial dependence of correlations was 
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observed out to a lag of about 0.5 m. Beyond 
this, correlations of grass species with SIR-C 
remained approximately constant. 


DISCUSSION 


Geostatistics can be applied to resource 
island data to provide several useful diagnostic 
features prior to actual mapping of the land- 
scape itself. For example, variography can 
define the presence and extent of spatial cor- 
relation and alert the researcher to apply with 
caution classical statistical comparisons that 
assume samples are independent and from 
identically distributed populations. For these 
methods to be more properly applied to spa- 
tial data, comparisons should probably be lim- 
ited to those samples separated by distances 
> range of the correlogram (Table 1; Webster 
1985, Robertson 1987). This is true for studies 
that compare samples collected along a con- 
tinuum such as distance, depth, or concentra- 
tion (i.e., resource gradient) or as a function of 
time. 

Another promising use of variography is to 
relate spatial continuity of two or more variables 
at the same site or the same variable at two or 
more sites by comparing variograms, covario- 
grams, or correlograms with one another. This 
approach may be useful for comparing the 
scale of ecological processes or ecosystem 
boundaries, but should be approached with 
caution for several reasons. First, each point in 
a traditional variogram represents the average 
value of the squared difference between many 
pairs of data points. While an average value 
may be appropriate for modeling spatial conti- 
nuity (as a summary statistic), it does not indi- 
cate the range of individual squared differ- 
ences or provide an estimate of the “goodness 
of fit” about each point in a variogram. The 
range of the deviation about the average value 
may be large or small (see Webster and Oliver 
1992), complicating comparisons of variograms. 
Thus, for comparative purposes other more 
“robust” representations of spatial dependence 
such as Journel’s mAD estimator, Cressie- 
Ilawkins’ robust estimator, or the rodogram 
(see Rossi et al. 1992) may be more appropriate 
choices. Second, even if variograms for two 
properties are sinilar, resultant estimates may 
vield very different maps (e.g., TI-N and SIR-C 
in this study). This is because estimation of un- 
known data values by kriging depends not only 


H20—C TI-N Vegetation 


SIR—C 
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Fig. 6.2 X 2-m maps of vegetation and kriged estimates of TI-N, H,O-C, and SIR-C for five ALE plots. Each kriged 
plot is composed of 1681 points estimated by ordinary point kriging. Vegetation maps indicate vertical projections of A. 
tridentata (black) and grass species (crosshatch) as determined from photographs. Soil properties are depicted in mg / kg 


soil (dw). 


upon a model of spatial continuity, but ultimate- 
ly upon degree and configuration of the known 
sample values in the field. 

During our analysis of resource island data 
using geostatistical methods, we made several 
assumptions or decisions about the data that 
could have affected our interpretations. First, 
we assumed that resource islands under A. tri- 
dentata could be monitored using a particular 
configuration of samples located within a 2 xX 
2-m plot. A different number of samples col- 
lected from a larger plot with a different shape 
or in a different pattern might have generated 
different correlogram models or kriged esti- 
mates (Webster and Oliver 1992). Second, we 
collected data at a single time during the year, 
thereby making the kriged maps “snapshots” 
in time and space. Data values and spatial 
continuity undoubtedly vary for some types of 
environmental variables (e.g., soil moisture or 


TI-N) on a seasonal or shorter time scale. 
Other environmental variables such as soil 
texture, total-N, or C might change more slow- 
ly. Third, we chose to analyze data for the five 
plots collectively rather than for each plot 
individually. This choice reflects an interpreta- 
tion that correlograms for individual plots were 
reasonably similar to each other and allowed 
our calculations to be based on a greater num- 
ber of paired comparisons. We reasoned that 
correlogram models of spatial continuity, de- 
rived from concatenated data, would stminarize 
typical patterns of spatial continuity and direc- 
tional anisotropies. Alternatively, although sin- 
gle-plot analyses would result in plot-specific 
models of spatial continuity with greater 
specificity, they might make generalizations 
difficult. Finally, we assumed that spatial con- 
tinuity for H»O-C was reasonably described 
by a single isotropic model, but we concluded 


39.4 GREAT BASIN NATURALIST [Volume 54 










A with Artemisia P B with grass species 
0.6 0 0 Oo 180° 
© 455, = 295. 
) 
0.4 Vv 90 Ete 270. 
< v135 4 315 
a 0.2 
° 
r 
“ 0.0 
"i _ 
° : 
: Saag 
-0.2 ised 
ia + 
oe 
-0.4 
S) 
| 
ro) 
N 
= 
re 
Q 
1 
“a 
"i 
° 
hee 
Oo 
oO 
| 
= 
n 
re 
Q 
1 
" 
"i 
o 
hes 
Oo 





0.00 0.25 0.50 0.75 1.00 1.25 0.00 0.25 Q.50 0.75 1.00 1.25 


Fig. 7. Directional cross-correlograms showing spatial covariance of indicator transformed TI-N, H,O-C, SIR-C data, 
and presence/absence (coded as 1 or 0) of A. tridentata and grass species. Correlograms were calculated for 0° (east), 
45° (northeast), 90° (north), 135° (northwest), 180° (west), 225° (southwest), 270° (south), and 315° (southeast) directions 
+30° tolerance. Each point shown summarizes a minimum of 5-4 pairs (range 54-158). 
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that directional anisotropies observed for TI-N 
and SIR-C were important enough to be 
accounted for in the estimation process. 

Results of a geostatistical analysis cannot 
completely replace “sound ecological reason- 
ing or theory (Rossi et al. 1992). Thus, the re- 
searcher must decide whether directional ani- 
sotropies observed in descriptive variography 
portray significant spatial patterns or are 
merely a coincidental result of the number 
and arrangement of data. The decision to 
account for spatial anisotropy in the kriging 
procedure is, in part, related to the desired end 
product of the geostatistical analvsis. For ex- 
ample, if the goal of an analysis is the most accu- 
rate representation possible of a particular 
resource island under a specific A. tridentata, 
then a highly detailed model of spatial continu- 
ity would be appropriate regardless of the 
source of spatial variability. In this case variog- 
raphy based on a concatenated data set might 
be less appropriate than analysis based only 
on the single plot. However, the goal of geo- 
statistical interpretation of ecological data may 
not be to produce detailed site maps. Instead, 
the ecologist may be more interested in pat- 
terns that are broadly applicable. Anisotropies 
can often be related to information about the 
environment such as stratigraphic, meteoro- 
logical, or hydrogeological patterns (Isaaks and 
Srivastava 1989) and may suggest linkages be- 
tween environmental variables. Our decision to 
account for anisotropies in the kriging process 
was, in part, influenced by information about 
another environmental parameter, prevailing 
wind direction. For this reason we would expect 
the anisotropies observed for TI-N and SIR-C 
to be a consistent feature of the ALE landscape. 

Directional correlograms revealed greatest 
spatial continuity for samples of TI-N and 
SIR-C in the 45° direction, northeast-south- 
west. Cross-correlograms, more specifically, 
indicated that above-median concentrations of 
soil properties were most correlated to A. tri- 
dentata in the 45° direction, or northeast. For 
the ALE site, cumulative records indicate that 
prevailing local wind direction is from the 
southwest quadrant (Table 3), corresponding 
to the downwind direction of greatest spatial 
continuity. Prevailing wind direction might 
influence spatial patterns of soil resources by 
affecting distribution of litter deposition which, 
for A. tridentata at the ALE site, may exceed 
60 kg / ha annuallv (Mack 1971). 
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TABLE 3. Frequency of eccurrence of wind at the ALE 
site. Source: H. Bolton, Pacific Northwest Laboratory. 


Source quadrant 1990 1991 
Hours Percent Hours Percent 
0-90° (NE) 1161 13.25 1208 13.97 
90-1S0° (SE) 1557 17 fo 1613 18.66 
1S0-270° (SW) 4090 46.70 3974 45.96 
270-360° (NW) aoe 2a 1851 21.41 


Evidence for the occurrence of resource 
islands in the ALE landscape was provided by 
comparing concentrations of soil resources 
collected near A. tridentata vegetation to 
those collected away from the plant. However, 
the specific sampling regime employed to 
evaluate “near” vs. “away influenced the par- 
ticular conclusion reached. For example, 
Bolton et al. (1990) were unable to conclude 
that concentrations of TI-N in soil under A. 
tridentata were significantly greater than con- 
centrations measured in open soil crust based 
on six samples drawn at random from each soil 
type (see also Doescher et al. 1984). In con- 
trast, we found that evaluating TI-N vs. dis- 
tance away from the A. tridentata axis resulted 
in the naive conclusion that significantly high- 
er concentrations of TI-N would always occur 
under A. tridentata plants (Fig. 3A, Table 1). 
Such a conclusion for TI-N and other soil 
properties would lead to a model of a land- 
scape composed of identically sized, symmet- 
rical resource islands centered on each A. tri- 
dentata individual and would infer some sort 
of causal relationship between concentration 
of TI-N and A. fridentata presence. However, 
kriged maps suggest that greatest concentra- 
tions of TI-N were not always associated with 
A. tridentata. 

Autocorrelation of soil properties was 
described using variography. The association 
of soil variables with A. tridentata individuals 
was supported jointly by kriged maps and 
cross-correlograms, with the latter showing 
that soil properties (especially Hj,0-C) were 
positively correlated to A. tridentata. Kriging 
is a means for producing visually satisfying 
maps of soil properties and provided addition- 
al insight into characteristics of resource dis- 
tribution under A. tridentata. However, we 
relied on these maps primarily as heuristic 
tools because we recognized that kriged maps 
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are models that can be influenced by decisions 
about the data set (e.g., concatenated vs. single 
plot), the “art” of variogram modeling, the 
type of kriging chosen (ordinary kriging is a 
data “smoother ), the specific search strategy 
used, and the method of graphical representa- 
tion. Finally, kriging by itself does not provide 
a measure of estimate confidence or reliability 
like nonparametric methods (Journel 1983) or 
stochastic conditional simulation (Rossi et al. 
1993). 

kKriged maps of H,O-C appeared to be 
most similar to graphs of summary statistics vs. 
distance from plant axis (Fig. 3A). In each kriged 
map (Fig. 6), high concentrations coincided 
with A. tridentata in a classic resource island 
pattern. These accumulations might be tied 
closely to inputs from A. tridentata litter fall 
representing a source of C that could be 
accessed by heterotrophic soil microorgan- 
isms. Alternatively, high concentrations of 
H1,0-C under A. tridentata might not indicate 
large C inputs. Instead, they might indicate 
the accretion of soluble, but reealcitrant, 
forms of C not readily useable by soil micro- 
organisms. In this case the term “resource 
island” would be ambiguous. To have ecologi- 
eal significance, a resource island must be 
evaluated for resource quantity, resource qual- 
ity, and presence of alternative resource sub- 
stitutes. Further, the significance of resource 
accuinulation into islands might change with 
time in relation to diurnal eycles, growing sea- 
son, or successional stage (Halvorson et al. 
1991). 

Kriged maps of SIR-C showed accumula- 
tions of soil microbial biomass in close proxim- 
ity to each A. tridentata individual. However, 
high concentrations were also observed for 
locations corresponding to other plant species, 
demonstrating that resource islands of micro- 
bial populations or aetivity can be numerous 
and are nonspecific to A. tridentata. Addition- 
ally, a significant amount of SIR-C was esti- 
mated for locations not associated with any 
plant. This suggests that while local inputs by 
plants may stimulate microbial population 
growth or activity, sufficient resources exist in 
the environment to support moderate amounts 
of SIR-C during some times of the year. 
However, plant loeation may control the dis- 
tribution of SIR-C indirectly through influ- 
ence on microclimatological factors such as 
soil temperature and evapotranspiration. 
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These factors would become more important 
during the bot, dry summer months and could 
limit distribution of SIR-C to locales closer to 
A. tridentata. 

Assessing the distribution of soil microbial 
populations or microbially mediated nutrient- 
cycling processes such as mincralization or 
denitrification is complicated by multiple 
resouree requirements and compensatory 
capabilities of living microorganisms (Smith et 
al. 1985). For example, microbial population 
size or activity within a C-substrate resource 
island might be limited by the availability of 
soil N. Conversely, the same microbial popula- 
tion might be limited by the availability of C- 
substrate despite an N-rich environment. 
Under such a scenario the greatest population 
size or activity might oceur in a location with 
low or intermediate quantities of both C and 
N, and the resource island for soil microorgan- 
isms or mineralization potential would appear 
distinct spatially from other resources. 
Estimation of soil properties like SIR-C that 
depend on the distribution of one or more 
other resources may need to be evaluated with 
respect to temporal and spatial distributions of 
alternative resources, 

Kriged maps of various soil properties can 
be interpreted within the context of the rela- 
tionship between the particular soil parameter 
and A. tridentata. Our data indicate that shape 
and orientation of resource islands under A. 
tridentata vary with the specific soil property 
considered, need not be centered on the axis 
of an A. tridentata plant, and need not be sym- 
metrical. The maps also provide evidence that 
suggests a vertical projection of the plant 
canopy is not well correlated to the distribu- 
tion of soil variables and thus should not be 
used as a basis for sampling designs (Fig. 6). 
For some soil properties (e.g., Hy0-C) the dif- 
ference between values characterizing the 
resource island and those characterizing the 
surrounding matrix may be large and the 
resource island may appear to have sharp 
boundaries. Conversely, the range of data val- 
ues for other soil properties (e.g., SIR-C) may 
be smaller and the transition from resource 
island to the surrounding soil matrix more 
complicated. Resource island boundaries may 
also change with direction, making sampling 
designs based on only a few transects ques- 
tionable. Finally, resource islands do not oecur 
under all A. tridentata or for all soil properties. 
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Other plants like annual and perennial grasses 
ean be the focal points for resource islands of 
some variables (Jackson and Caldwell 1993b). 

Geostatistics allows estimation and map- 
ping of resource islands in considerable detail. 
Such maps can be used to further our under- 
standing of the ecology of A. tridentata, refine 
nutrient budgets for shrub-steppe ecosystems, 
reveal the existenee of resource and process- 
dependent patterns, and help provide a ratio- 
nale for sampling designs based on natural 
boundaries. Besides two-dimensional space, 
geostatisties ean be used to consider differ- 
ences in spatial continuity with soil depth (i.e., 
a third dimension) or time (via repeated mea- 
surements). Hlowever, even with geostatisties, 
our definition of a resource island can be 
improved. Whether a resource island is more 
properly delineated by some minimum differ- 
ence in resouree coneentration or related to 
the eeologieal significance of small differences 
in concentration remains to be answered. 
Further, the resouree island “effeet” may be 
related to more than a single environmental 
parameter. Consequently, methods must be 
developed to simultaneously integrate infor- 
mation for several environmental variables 
and summarize them spatially. 
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